Single-neuron criticality optimizes analog dendritic computation 
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Neurons are thought of as the building blocks of excitable brain tissue. However, at the single 
neuron level, the neuronal membrane, the dendritic arbor and the axonal projections can also be con- 
sidered an extended active medium. Active dendritic branchlets enable the propagation of dendritic 
spikes, whose computational functions, despite several proposals, remain an open question. Here we 
propose a concrete function to the active channels in large dendritic trees. By using a probabilistic 
cellular automaton approach, we model the input-output response of large active dendritic arbors 
subjected to complex spatio-temporal inputs, and exhibiting non-stereotyped dendritic spikes. We 
find that, if dendritic spikes have a non-deterministic duration, the dendritic arbor can undergo 
a continuous phase transition from a quiescent to an active state, thereby exhibiting spontaneous 
and self-sustained localized activity as suggested by experiments. Analogously to the critical brain 
hypothesis, which states that neuronal networks self-organize near a phase transition to take ad- 
vantage of specific properties of the critical state, here we propose that neurons with large dendritic 
arbors optimize their capacity to distinguish incoming stimuli at the critical state. We suggest that 
"computation at the edge of a phase transition" is more compatible with the view that dendritic 
arbors perform an analog and dynamical rather than a symbolic and digital dendritic computation. 



Critical systems are organized in a fractal-like pattern 
spanning across numerous temporal and spatial scales. 
The brain has recently been included amongst abundant 
physical and biological systems exhibiting traces of crit- 
icality [l|, y. In the past decade, several phenomena 
suggestive of critical states have been observed in dif- 
ferent systems: multielectrode data from cortical slices 
in vitro |3H5|, anesthetized [fj, awake [7[ and behaving 
animals [8[; human electrocorticography |9|, electroen- 
cephalography [lOj, magnetoencephalography IIOL llll l , as 
well as functional magnetic resonance imaging [111 . Il2j 
recordings. Together, these experiments give robust ev- 
idence for the intrinsic pervasiveness of criticality into a 
wide range of spatio-temporal brain scales. 

The critical brain hypothesis offers an appealing solu- 
tion to a long-lasting conundrum of neuroscience: how 
can localized information (sometimes from very spe- 
cific regions) propagate in the brain without a spa- 
tial/temporal exponential decay or saturation [l|, Il3l |? 
Advantages of neuronal systems poised around a crit- 
ical state include optimization of the dynamic range 
of neuronal networks 1J|, as confirmed experimentally 



[15j . as well as transmission and storage of informa- 
tion 0, |j, LLs, LL7[ . However, a fundamental assumption 
of the critical brain hypothesis has so far not been ex- 
amined: is the single neuron the minimal dynamic unit 
in neuronal networks, like a spin site in the Ising model 
for ferromagnets (the prototypical system for phase tran- 
sitions)? Or can critical phenomena associated with 
phase transitions already occur at the single neuron level? 
Adopting a statistical physics approach, we explore the 



behavior of dendritic branchlets at the sub-cellular level 
as forming a network of functional dynamical units. 

Here we treat neurons with their extensive dendritic 
trees as excitable media (see Fig. QJi). This task is far 
from trivial because of the complexity [18j, |l9j, the di- 
versity of neurons [20|, and the absence of information 
about key elements governing the dynamics, such as how 
the zoo of ionic channels is distributed along the den- 
drites |2ll424| . In fact, there is an entire research area 
concerned with biophysically detailed neuronal model- 
ing [24H27 1 and dendritic computation 28 1 . 

For our modeling purposes, we evoke the statistical 
physics principle that emergent phenomena seem to de- 
pend on only few characteristics such as symmetries, 
dimensionality, network topology, type of coupling etc. 
Hence the basic dynamical units (atoms, spins) need not 
be modeled in detail [29(. Our cellular automaton ap- 
proach to dendritic computation unveils the emergence 
of critical phenomena at the single neuron level from a 
large number of sub-cellular interacting units (dendritic 
compartments) . 

We find that the dendritic arbor can fire spontaneously 
when we take into account the fact that dendritic spikes 
are non-stereotyped, that is, with variable durations [30(. 
In this case, a smooth continuous phase transition ap- 
pears between rest and self-sustained dendritic activity. 
At the interface between these states, there is a critical 
regime, which optimizes the dynamic range of the input- 
output response function of the neuron's dendritic arbor. 
The presence of the phase transition can increase the 
maximum neuronal signal compression capacity by up to 
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FIG. 1. Model of an active dendritic tree with non- 
stereotyped dendritic spikes, (a), Excitable elements 
(circles) connected (bars) in a Cayley tree topology with 
G = 3 layers and coordination number z = 3 (one mother and 
k = 2 daughter branches), (b)-(d), Biologically motivated 
dendritic spike with duration depending on the distance from 
the soma as a net effect of variable density of ionic channels. 
Red arrows indicate that the site is being stimulated (by an 
external input and/or by mother or daughter branches), and 
s represents the average period of dendritic spikes at layer 



~ 100 times. In our critical-neuron model, the mini- 
mal dynamic units are dendritic branchlets (or perhaps 
patches of dendritic spines), so that even a single neu- 
ron can have a highly nonlinear input-output response 
function. This capability of a single neuron to compress 
stimulus intensity varying over several orders of magni- 
tude in a decade of output firing rate could be the basis 
of psychophysical power laws [3l[ . In addition, such re- 
finement in the basic dynamical unit corresponds to a 
modeling improvement of spatial resolution of a few or- 
ders of magnitude. 



RESULTS 

The overwhelming majority of models of criticality in 
neuronal networks neglect the characteristic tree shape 
of neurons. Typically, the collective behavior studied is 
generated by oversimplified point-like neurons with no 

a fir 



spatial structure [4j, [14|, I32H361]. Dropping this strong 
assumption unveils a new sublevel of neuronal dynamics 
where concepts from statistical physics can be applied. 
In previous work we introduced a simple model of active 
dendrites [371] . In the model, nonlinear excitable waves of 
activity, not accounted for by passive dendrites and linear 
cable theory, propagate and annihilate upon collisions, 
giving rise to large signal compression abilities. Weak 
inputs are highly amplified whereas strong inputs are not 
subjected to early saturation. The dynamic range, which 
quantifies how many decades of stimulus intensity can 
be distinguished by the excitable medium, attains large 



values. This property has proven robust against several 
model variants [371.138 1. 

In fact, excitable media with a tree structure per- 
formed better than other network topologies [39|, |40J. 
Moreover, the system dynamics in a tree topology al- 
lows analytical treatment due to the absence of loops. 
The model was quantitatively understood by solving the 
master equation under an excitable-wave (EW) mean- 
field approximation [38(. However, phase transitions are 
not observed in such dendritic arbors when the dynamics 
of the excitable elements (active dendritic branchlets) is 
deterministic. In this paper, we study a different den- 
dritic arbor with non-stereotyped dendritic spikes that 
exhibit probabilistic and non-homogeneous (layer depen- 
dent) spike duration. Such biologically motivated im- 
provements in the model give rise to a new dynamical 
regime consisting of a self-sustained state with sponta- 
neous activity that emerges through a non-equilibrium 
phase transition. 

Physiological experiments show that dendritic spikes 
depend on several voltage-gated ionic channels. Due 
to the non-homogeneous density of channels along the 
proximal-distal axis, dendritic spikes can present distinct 
dynamical features in different regions of the dendritic 
tree. Dominated by the slow calcium-dependent poten- 
tials [27IJ41I E2] or N-methyl-D-aspartatc (NMDA) po- 
tentials [30j , the dendritic spikes/plateaus become longer 
at more distal sites, as shown in Figs. [Tp through d. 
These responses account for a variety of active behaviors 
including dendritic spikes with longer and more variable 
duration than sodium-dependent somatic spikes. 

In our previous models 37], |38[ , other parameters have 
been studied, but the time- length of the spikes that we in- 
troduce here proves to be a crucial factor. Such variabil- 
ity in the duration of the spike has been shown to shape 
the network dynamics [43j , and to enhance the computa- 
tional capability of neuronal networks }44| . In the present 
model, the variable duration of dendritic spikes plays 
an important role in dendritic computation because it 
is capable of controlling the emergence of a phase tran- 
sition. Evidence for both sides of this continuous non- 
equilibrium phase transition has been reported experi- 
mentally [23j, |41[ . Among the active states there are sev 
eral different pat terns of activity: plateaus [3 
spikes [HI 0,113, bursts [H,|l9J, and oscillations 
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As illustrated in Fig. QJi, we consider the dendritic ar- 
bor as a probabilistic cellular automata with a Cayley 
tree topology of active branchlets capable of transmit- 
ting dendritic spikes. The states update in parallel with 
discrete time steps of St = 1 ms. Each excitable unit 
(the i-th dendritic branchlet) has three possible discrete 
states Si(t) £ {0,1,2}. A quiescent site (si(t) = 0) may 
become active in the next time step (si(t+ 1) = 1) either 
by external driving or by propagation of an excitable- 
wave from an active neighbor with probability p\. Ex- 
ternal input arrives at quiescent sites with probability 



Ph = 1 — cxp(— hSt) per time step, where h stands for 
the rate of independent Poisson processes and may vary 
in a range of several orders of magnitude. An active site 
(si(t) = 1) becomes refractory (sj(i + 1) = 2) with prob- 
ability pg at each time step, where g indicates the layer 
that the site belongs to (see Fig. [Tb) . To close the cycle, 
refractory sites become quiescent with probability p 1 (we 
fix p 1 = 1/2 throughout this paper). 

The model generalizes our previous work |37j, l2S| by 
introducing the probabilistic spike duration. In what 
follows, we first consider a simple homogeneous model 
with spike duration controlled by a uniform probabil- 
ity Ps — Ps < 1, V<? (our previous model [37J is thus 
recovered when ps = 1). For ps ~ the spike dura- 
tion becomes unreasonably long. Thus, our attention is 
mostly concentrated in the regime between the extremes 
(0 < ps < 1). This homogeneous case has some advan- 
tages: it is clear, intuitive, and allows analytical insights 
about its phase transition, which is controlled by p\ and 
Ps- 

Next, we study a more realistic scenario based on the 
relevant physiological evidence previously described: the 
model assumes that the spike duration depends on the 
apical distance. We consider a linear dependence p$ — 
1 — 0. 9-jka, where G = 10 (unless otherwise stated) stands 
for the tree size. The parameter a (0 < a < 1) controls 
the inhomogeneity, and it is rather strong for a < 1. 

To characterize the dynamical regimes of our system, 
we must define our order parameters [29j. The main or- 
der parameter of interest here is the stationary firing rate 
of the apical site so (see Fig. [T^,): 



is reached. If the tree reaches an active and stable level, 
the system is said to be supercritical. On the contrary, 
if the activity of the system vanishes rapidly, the system 
is subcritical. For large systems, the fate of the system 
depends weakly on the initial condition and is mainly 
determined by control parameters p\ (which controls the 
coupling between branchlets) and ps (which controls the 
average spike duration). The critical regime occurs at 
the border between the sub- and supercritical regimes, 
where the activity level vanishes slowly and without a 
characteristic time scale [14|. 
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where T is an averaging time window and 5ij is the Kro- 
neckcr delta (= 1 if i = j, zero otherwise). The apical 
dendrite firing rate F is biologically interesting because 
it provides the main input to the neuronal soma and can 
be thought of as proportional to the neuronal firing rate. 
Notice that so and hence F are functions of the external 
input with intensity h. Our main focus will be on the 
response function F(h). 



Phase diagram without external input 

It is instructive to analyze the complex spatio-temporal 
dynamics of an extensive dendritic tree in the simplest 
scenario, when average spike durations are homogeneous 
(a = 0) and the system is subjected to no external driv- 
ing (h = 0). The absence of stimulus is a somewhat artifi- 
cial condition, but very informative because it crisply un- 
covers the critical region. Starting with a small fraction 
of active nodes, we follow the dynamics for a sufficiently 
long time (T ~ 10 4 time steps) until a stationary state 



FIG. 2. Continuous phase transition in active dendritic 
trees with non-stereotyped dendritic spikes, (a), Av- 
erage firing rate as a function of p\ (which governs the cou- 
pling among branchlets) and ps (which governs the average 
spike duration) in numerical simulations for tree size G = 10. 
(b) , Returning probability (see text for details) as a function 
of p\ and pg. (c), Same as panel (a) but for the general- 
ized excitable-wave mean-field approximation (see Methods 
for details). 



As represented in Fig. [5^, we numerically find the crit- 
ical curve of the parameter plane (p\,ps) in a finite tree 
with G = 10. The critical curve corresponds to the bor- 
der between the absorbing state with F(h = 0) = (blue) 
and the active state with F(h = 0) > (red). Naturally, 
a true phase transition occurs only for infinite systems, 
but for moderately large systems we already observe ac- 
tivity that persists for simulation times much longer than 
any relevant biological time scale. It is important to em- 
phasize that no active phase exists for dendritic spikes 
with a deterministic duration, p$ = 1. Only when they 
are non-stereotyped (ps < 1) does self-sustained activity 
become possible. 

To qualitatively understand why the activity dies out 



in a specific part of the parameter space of a finite tree, 
we define the returning probability R. It corresponds to 
the probability that an active site A stimulates a given 
quiescent neighbor B and receives back the stimulus at 
later times after performing a complete cycle (i.e., after 
going through the refractory (sa = 2) and the quies- 
cent (sa = 0) states). A necessary (but not sufficient) 
condition for a network without loops to exhibit a stable 
self-sustained state (and consequently a phase transition) 
is that R =/= 0. A nonzero returning probability is neces- 
sary for the persistence of the self-sustained state because 
otherwise the activity ceases after a few time steps due 
to collisions of the excitable waves with the boundaries 
(layer G) and with one another. The Methods sectiort il- 
lustrates the fastest example of a successful process of re- 
turning activity, which requires at least three time steps, 
and extends the calculation for arbitrarily long two-site 
processes. Figure [2b shows R(p\,ps), in which a transi- 
tion similar to that of Fig. [2k is observed. A comparison 
between them confirms that R ^ is a necessary but not 
sufficient condition for an active phase to emerge. 

Another analytical approach to shed light on the re- 
sults is to employ one of the many mean-field approxi- 
mations available in the statistical physics literature [29| . 
We have previously developed an excitable-wave (EW) 
mean-field approximation that focuses on the wave prop- 
agation direction, making it particularly suitable for ex- 
citable waves on a tree [38|. In the Methods section we 
describe a generalized version of the approximation to 
account for non-stereotyped dendritic spikes. As Fig. [2fc 
shows, the phase transition observed in the simulations 
can be qualitatively described by the generalized EW 
(GEW) mean-field approximation. 



Neuronal response to external input 

From herein, we focus on the more relevant scenario in 
which a neuron has to somehow cope with information 
arriving stochastically at its thousands of synapses. We 
make the simplest assumption that, at each branchlet, 
the average excess of incoming excitatory post-synaptic 
potentials (EPSPs), as compared to inhibitory post- 
synaptic potentials (IPSPs), can be modeled by an in- 
dependent Poisson process with rate h [52j. We study 
the response function F(h) (averaged over T = 10 4 time 
steps and five realizations) and its dependence on model 
parameters. 

Figure^ depicts three response functions F(h;p\) for 
a fixed ps, exemplifying one response function of each 
kind: subcritical (triangles, blue), critical (closed cir- 
cles, black) and supercritical (open circles, red). Al- 
ternatively, for a fixed p\, the critical line can also be 
crossed by varying ps, as depicted by Fig. [3b- In this 
case, the maximum firing rate displays its dependence 
on p s : F max = (1 + ips)^ 1 ■ 
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FIG. 3. Response functions and phase diagram, (a), 

Response curves F(h) for p\ = 0.5, 0.755, 1 (from right to 
left). Horizontal and vertical arrows show the relevant pa- 
rameters for calculating the dynamic range A (see Eq. ©). 
(b), Response curves dependence on ps for homogeneous dis- 
tribution p 3 s = pg^g. Symbols connected by dashed lines 
represent simulation results whereas continuous lines repre- 
sent the generalized excitable-wave (GEW) mean-field ap- 
proximation, (c) and (d), Response curves dependence on 
g, i.e., p 3 — 1 — 0.9-j^a. (d), Family of response curves for 
Px = 0.4, 0.6, 0.8, 0.85, 0.9, 1 (from right to left). Inset: phase 
diagram. Tree size is G = 10 for all panels. 



For some classes of dendrites, it can be more realistic 
to employ our heterogeneous model, with a > 0. We con- 
sider a varying in a logarithmic range of several decades 
because the phase transition already occurs for a <c 1. In 
this instance, the phase transition is controlled by both 
p\ and a, giving rise to a critical curve in parameter 
space, as depicted in the inset of Fig. (3Ji. Analogously 
to the homogeneous case, for fixed p\ the transition de- 
pends on a, as depicted in Fig. [3J:. For both homoge- 
neous and heterogeneous models (solid line of Figs. [3b 
and c) , the GEW mean- field approximation captures par- 
ticularly well the behavior of both sub- and supercritical 
curves but fails to capture the strong amplification of the 
critical curves under weak external driving. As shown 
in the family of curves for fixed a displayed in Fig. [3ji, 
the critical curve has a very small exponent (m ~ 0.11), 
which cannot be correctly described by mean-field ap- 
proximations. This exponent m is similar to a Stevens 
psychophysical response exponent [311 ] observed at the 
single- neuron scale [37[, as had been previously noticed 
at the level of a neuronal ensemble [1J, [53J, [54[ . 

Excitable networks have a recognized ability to com- 
press several decades of input rate in a single decade 



of output rate [lj, |36Tt4Q| . [53j. This information pro- 



cessing capacity emerges exclusively from local inter- 
actions between the excitable elements. A large dy- 
namic range is robust, being also obtained for mod- 
els with increasing levels of biophysical realism, such as 



networks of FitzHugh-Nagumo and Hodgkin- Huxley ele- 
ments 53j,|55|, as well as more sophisticated models based 



56,57 



on detailed anatomical information of the retina 

The definition of dynamic range is illustrated in 
Fig. [5k. [14[. We start by neglecting the regions of the re- 
sponse curves which are close to the detection threshold 
F m in or the saturation level F max . According to a con- 
ventional definition, only the interval between Fo.i an d 
F0.9 (arbitrarily defined within 10% of the plateau levels) 
properly codes the input. Those firing rates correspond 
to external driving with rates ho.i and ho.g. The dy- 
namic range is the range (measured in decibels) between 
/lo.i and /iq.9j i-e., 



v «0.1 



(2) 



The phase transition is an important feature for signal 
compression. Comparing the different regimes, as de- 
picted for instance in Fig. [3^, a simple reasoning allows 
one to understand why the critical curve optimizes the 
compressing capacity 14[ . The subcritical curves cannot 
amplify small external stimuli, since the activity from in- 
coming pulses vanishes exponentially fast in space and 
time. At the other extreme, for supercritical curves, al- 
most any pulse initiates the self-sustained network ac- 
tivity for arbitrarily long time periods, masking the re- 
sponse to weak stimuli. However, when the network has 
the correct critical parameters, the spontaneous network 
activity is composed of neuronal avalanches launched by 
spontaneous fluctuations of the system. In the presence 
of an external field h, such avalanches add and super- 
pose, creating the response function F(h) (analogous to 
magnetization for non-zero fields in magnetic systems). 
In this case, the network dynamic range is optimized. 

Since phase transitions are properly defined only for in- 
finite systems, we have investigated the effects of system 
size in our model. It turns out that, even for moder- 
ately large systems (e.g. G = 10, as in Figs. [2] and [3]), 
self-sustained activity survives for a period of time which 
is much longer than any relevant characteristic time for 
neuronal processing (T ~ 10 4 time steps ~10 seconds), 
which is the operational marker that we have employed 
for the transition. Its occurrence is associated with the 
peaks in dynamic range shown in Fig. 0^. 

In small trees, activity dies out rapidly due to wave 
collisions with the boundaries, preventing self-sustained 
activity. Taking for example the case of G = 4 in Fig. 2k, 
the maximum A max (of the A vs. p\ curve) occurs at 
px = 1, meaning that finite size effects masks the (infi- 
nite size) phase transition. Larger trees, in contrast, can 
exhibit A max = A(p\ c ) for weaker coupling revealing the 
presence of a phase transition at p\ c (G) < 1. The larger 
is the tree, the smaller is the p\ c (G) (Fig. 0^), and the 
larger is the A max (Fig. gb). 

As presented in Fig.[5J pg = 1 prevents an active state. 
However, for p$ < 1 a phase transition to an active state 



occurs in large enough trees. The black arrows in Fig. [4]b 
represent the minimum tree size that gives rise to our 
operationally defined transition for ps = 0.5 (G = 4), 
&ndp s = 0.9 (G = 6). 
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FIG. 4. Dynamic range A. (a) and (b), for homogeneous 
distribution of pa. (c) and (d), for layer dependent distribu- 
tions: pg = 1 — 0. 9^a. (b) and (d) depicts A maa! which is the 
maximum of the A(p^) curve. In each curve phase transitions 
occur only at and on the right hand side of arrows. Sym- 
bols connected by dashed lines represent simulation results 
whereas continuous lines represent the generalized excitable- 
wave (GEW) mean-field approximation. 

The optimization of the dynamic range at the critical 
value (Fig. |4j:), and the growth of the maximum dynamic 
range with the tree size (Fig.|4ji) are also observed in the 
model with non-uniform spike duration. In this more re- 
alistic model version, the phase transition vanishes for 
a ~ 0, as depicted in the inset of Fig. [3ji, since in such 
case the spike duration is deterministic across the tree. 
Increasing a from zero, A mQa; (a) grows up to a maximum 
around the region where the phase transition appears 
(represented by arrows in Fig. gjl). As shown in Fig. gji, 
larger values of a lead to a phase transition with smaller 
critical coupling p\ c . On the other hand, as depicted in 
Fig. HJi, the A„ iax (a) curves show plateaus with heights 
and widths that both increase with tree size. This sug- 
gests that large active dendritic trees with non-uniform 
spike duration lead to a large dynamic range; a result 
whose robustness is attested by the width of the plateaus. 
In particular, the difference in dynamic range for an op- 
timized a compared with the a = case is remarkable, 
attaining about 20 dB for G = 5 (Fig. gH). 



DISCUSSION 

We have shown that the idea of a critical network of 
excitable elements being able to optimally process incom- 
ing stimuli can be applied at the subcellular (dendritic) 



scale. Instead of a network of excitable point-like neu- 
rons connected by synapses, here we have shown that 
electrically connected excitable dendritic branchlets can 
cause a dendritic arbor of a single neuron to undergo a 
phase transition. The key ingredients for this scenario 
are stochasticity in the duration of dendritic spikes and 
stochasticity in the propagation of the excitation from 
one branchlet to its neighbors. 

Several experiments indicate that a dendritic spike at 
a distal site is usually not strong enough to trigger a 
somatic action potential [58(. This suggests that real 
dendrites ought not typically show deterministic signal 
propagation (that is, p\ < 1). In the present model, tak- 
ing for instance a dendrite with tree size G = 10 and 
probability of activity propagation p\ — 0.5, the proba- 
bility for a most distal dendritic spike to propagate all 
over to the apical site and potentially generate a spike 
is very small: 0.5 10 ~ 10 -3 . Thus, the model is in ac- 
cordance with experimental wisdom. Remarkably, every- 
when the dendritic parameters lie near the silent/active 
phase transition, our results prove that non-reliable den- 
dritic spike propagation is compatible with an optimal 
dynamic range. 

We have examined the implications of non-stereotyped 
(non-homogeneous and non-deterministic) dendritic 
spikes for dendritic computation. We have shown that 
a necessary condition for a phase transition to occur 
in active dendrites is that the returning probability R 
be nonzero, which in turn can only occur with non- 
stereotyped dendritic dynamics (ps < 1). In an attempt 
to incorporate more details on dendritic spikes, we have 
also modeled different spatial functions of average spike 
duration, and showed that our main results are robust in 
that parameter space. 

While the non-stereotyped criterion must be satisfied 
to give rise to a critical neuron, it docs not play such 
a crucial role in giving rise to a critical neuronal net- 
work [431] . The key difference between the two spatial 
scales is the topology: due to inhibitory signaling mech- 
anisms during growth [59[, a dendritic tree contains no 
loops, whereas in neuronal networks they abound. The 
possibility of criticality at both scales naturally raises 
the question: if neurons can be critical, why should we 
need critical networks? We propose that the answer de- 
pends on the type of neuronal computation and the scale 
at which it occurs. In the cases where the processing is 
distributed among a large number of neurons with small 
dendritic trees, it is plausible to look at point-like neurons 
as the basic units of a larger system which, via balanced 
synapses, can tune itself to a collective critical regime. By 
contrast, single- neuron criticality could play a computa- 
tional role in cases where neurons have extended den- 
dritic trees and therefore must cope with large variations 
of synaptic input (such as olfactory mitral cells or cere- 
bellar Purkinje cells, for instance). Evidently, in princi- 
ple nothing prevents the mechanisms at both scales from 



acting together. 

The input-output response function of a critical neuron 
amplifies small-intensity inputs while preventing early 
saturation when subjected to large-intensity inputs. The 
response function follows a slow-increasing power-law 
function F ~ h m , in which m corresponds to a rather tiny 
exponent m ~ 0.11, as shown in Fig. [3ji. Such a slow- 
increasing function could be easily confounded with a log- 
arithmic (Weber- Fcchncr's psychophysical) law 3l[ |53|. 
These results strengthen previous suggestions that a sys- 
tem poised at a critical point is a natural candidate 
to explain how m < 1 Stevens' power-law exponents 
emerge 14 1. 

Our simple model yields similar results to those ob- 
tained with detailed compartmental modeling of mor- 
phologically reconstructed dendrites: the maximum dy- 
namic range a dendritic arbor can attain grows with the 
tree size [57] . This supports the previous proposal that 
neurons with large dendritic arbors might have grown 
to attain such impressive size and complexity in order 
to enhance their capacity to distinguish the amount of 
external driving 37J, |38| . 

Owing to our irresistible tendency of comparing brain 
processes with whichever technology happens to be dom- 
inant at the time, the term dendritic computation is typ- 
ically associated with a symbolic and digital-like infor- 
mation processing [2l|, |60(. However, dendritic trees are 
living tissues that change all the time, growing and re- 
tracting branchlets, spines and synapses. For example, 
30% of spine surface retracts in hippocampal neurons 
over the rat estrous cycle [6l(. This, we believe, does not 
seem compatible with a fixed circuitry implementing, say, 
logical gates. In contrast, our framework suggests that 
dendritic arbors perform a robust analog computation, 
which is resilient to disturbances of tree properties. For 
instance, pruning half of the tree corresponds to chang- 
ing from G to G — 1 , which amounts to a small decrease 
in the dynamic range. 

In contrast to physical systems, biological systems can 
approach the critical state through homeostatic mech- 
anisms [33j, |62| . This self-organized tuning is essential 
for processing information over a large range of stimulus 
intensities. We believe that critical neurons do indeed 
exist, but experimental limitations or lack of a clear the- 
oretical framework possibly prevented spatio-temporal 
neuronal criticality from having been reported in the lit- 
erature to date. 

Experimental confirmation of single-neuron criticality 
requires spatial and temporal recording resolutions at the 
edge of current available techniques, which are likely to 
be available in the near future. Several lines of evidence 
already suggest critical phenomena at the neuron level. 
In vivo long-range inter-spike time interval correlations, 
which may have originated from a critical neuron, have 
been reported in human neurons from hippocampus [63j 
and amygdala [64j . as well as in rat neurons from the 



leech ganglia and the hippocampus 65|. Furthermore, 
neglecting the extensive spatial features of neurons, crit- 
icality in neuronal excitability has recently been pro- 
posed [351 ] - 

Signs of a self-sustained state can be associated to pat- 
terns of dendritic activation: spikes, plateaus, bursts, 
and oscillations [4l|. Evidently, dendritic trees can also 
be on the other side of the spectrum, showing a silent 
or rest state. Our model predicts a continuous phase 
transition so that criticality lies between such subcritical 
(silent) and supercritical (self-sustained active) regimes. 
We speculate that adaptation and homeostatic mecha- 
nisms, similar to those discussed in Refs. [33j, |62j, could 
finely tune a self-organized critical state (or more pre- 
cisely, a self-organized quasi-critical state [62j). Such 
homeostatic mechanisms, conceivably present at the den- 
dritic level, will be presented in a forthcoming paper. 

On the computational side, our cellular automaton 
model can be generalized so that each branchlet is 
modeled by a detailed biophysical compartment with a 
plethora of ion channels [57| , coupled with (noisy) axial 
resistances. However, as in the proverbial arbor that pre- 
vents us from seeing the forest, we must be aware that de- 
tailed biophysical modeling may hamper us from detect- 
ing phase transitions if the model has few compartments. 
Our statistical physics model, with simple excitable units 
but connected in a large dendritic network, enables us to 
see the forest. 



METHODS 
Returning probability calculation 

The returning probability R > is a necessary con- 
dition for a phase transition. For the phase transition 
between quiescent and self-sustained states to occur the 
system must allow the self-sustained activity to be kept 
for arbitrarily long times. In a finite network without 
loops, such as a Cayley tree, this condition is only satis- 
fied if the activity can go back and forth between two 
neighbor sites, say A and B, as illustrated in Fig. [5] 
Otherwise, in the absence of external driving, the sys- 
tem must go to the rest state (regardless of the initial 
condition) after a maximum of 2G + 1 time steps. 

It is possible to calculate the returning probability be- 
tween two neighbor sites that obey the cyclic cellular 
automata rules. The example of Fig. [5] can be com- 
puted directly. Considering the initial condition at time 
t, SA(t) = 1 and sb (t) = 0, there is a unique possible path 
for the excitable wave to go back and forth after three 
time steps. Notice that, by construction, SA(t) must com- 
plete the cycle (i.e., it must go through states sa = 2 
and sa = 0), since we have excluded the trivial solu- 
tion: SA(t) = 1, Vt. The final configuration s^(i + 3) = 1 
and Sb(£ + 3) = 2 is therefore reached with probability 
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FIG. 5. Example of returning activity. An active site 
A stimulates its neighbor B (red arrow) and receives back 
the activity three time steps later. This process involves a 
series of intermediate steps. Site A must pass through the 
refractory (sa(£ + 1) = 2) and quiescent (sa(£ + 2) = 0) states 
before becoming active (sa(£ + 3) = 1) again. Moreover, site 
B must be kept in the active state (ss(£ + 1) = 1) for at 
least one time step (ss(£ + 2) = 1) in order to be able to 
excite back the susceptible A site (sa(£ + 2) = 0). Therefore 
a nonzero persistence probability (1 — ps > 0) is fundamental 
for a nonzero returning probability (R > 0). 



p|(l — Ps)p\p^, where psP-yPx is the probability of site A 
to go throughout the cycle, and p\(l — ps)ps is the prob- 
ability of site B to receive the input, to remain active, 
and finally to become refractory. 

This example is the simplest one. In general, how- 
ever, there are infinitely many possibilities that must be 
included in a full calculation of R. The configuration 
,s_4 = 2 and sb = 1, which was depicted at time t + 1 
in Fig. [SJ could have also be found at later times: t + 2, 
t + 3 and so on, leading to identical final configuration 
sa = 1 and sb = 2. The calculation of the returning 
probability which accounts for all the possible intermedi- 
ate steps, considering that site B has just been activated 
(ss(t + 1) = 1), involves three sums of infinite geometric 
series, Si, S 2 and S3. Site A can spend arbitrarily long 
periods active (oc Si), refractory (oc S 2 ), or quiescent 
(oc S3) and nevertheless receive back the excitable wave, 
given that the pivot site B remains active in the mean- 
time. Finally, without loss of generality, we can ignore 
the final state of site B to find the returning probability 
after infinitely many time steps: 



where 



R = psp~/(l - ps)p 2 \SiS 2 S 3 

1 
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Equations ©-(JU) were used in Fig. QJ. However, this 
result can also be extended to account for the heteroge- 
neous case with ps belonging to different layers (a and 



b). The returning probability R a,b turns out to be 

b\ 2 rta,b ria,b rya,b 
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Generalized excitable-wave mean-field 
approximation 

This approximation generalizes the recently proposed 
excitable- wave (EW) mean-field calculation |38[. Here, 
the scope of the approximation is enlarged to account for 
non-stereotyped dendritic spikes with probabilistic spike 
duration (ps < 1) and non- homogeneous spatial distribu- 
tions (ps(g))- This is a single-site mean- field approxima- 
tion, but it keeps track of the excitable- wave direction of 
propagation. Remarkably, the system response is better 
captured by this approximation than by the traditional 
two-site mean-field approximation [38j |. 

First we explain the notation and then recall the sys- 
tem master equations. Assuming that at time t a site at 
generation g is in state x; its mother site at generation 
g — 1 is in state y; and i (j) of its daughter branches at 
generation g+1 are in state z (w) etc., the joint probabil- 
ity of this configuration reads as: Pf (j/; x; z' 1 ' , w^> ,...). 
We also employ the usual normalization conditions 

Y, W3 p t (y> x > w i' w 2' w 3) = p i ( x '>y> w ii w 2)- Thus ' tlic 

master equations for arbitrary coordination number z = 
k + 1 and layers < g < G is given by [38( : 



P/ +1 (;l;) = P/(;0;)-(l-p fe ) 



E 



](-irP/(;0;lW) 



i+ 



ifk 



(-lYPf l;0;l(' 
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+(l-p s )Pf(;l;), 
P t s +l (; 2; ) = PS P?(; 1; ) + (1 - p 7 )/f(; 2; ) , 
P/ +1 (;0;) = l-Pf +1 (;l;)-P/ +1 (;2;), 



(11) 
(12) 
(13) 



where Pf(y; x; w^ ') = Pf(y; x; ) is a two-site joint prob- 
ability. For simplicity, we also drop the excess of notation 
symbols (;), such that Pt(x) stands for P t 9 (; x; ) which is 
the probability of finding at time t a site at generation g 
in state x (regardless of its neighbors). 

Equations controlling the active state for sites belong- 
ing to the first (g = 0) and last (g = G) layer can be 



obtained from straightforward modifications of Eq. (|lll) , 
yielding: 

f? +1 (;l;) = i?(;0;)-(l-p h ) 

fe+lr 'fc + 1 



5>i [*y (-irp t °(;o ; i«) 



+(l- TC )P t °(;l;), (14) 

P t G +1 (; 1; ) - P t G (; 0; ) + (1 - p fc ) [/3p A P t G (l; 0; )] 

+(l-p s )P t G (;l;), (15) 

while equations controlling the refractory (|12[) and quies- 
cent (|13p states remain unchanged. The full description 
of the dynamics requires higher-order terms (infinitely 
many in the limit G —> oo), but, as we show below, 
Eqs. ([Ti"l) - (fT3")) are enough for the GEW approximation. 

g>0 




FIG. 6. Generalized excitable-wave mean-field approx- 
imation. Top- left panel illustrates the dynamics of each node 
and of layer g = under the GEW mean-field approximation. 
Top-right panel illustrates the dynamics of further layers un- 
der the GEW mean-field approximation. Bottom panel shows 
the excitable-waves direction of propagation. Notice that in 
the absence of purple arrows the flux vanishes in the state 
(1C) of the last layer. The purple arrows represent transi- 
tions that occur at each time step with probability (1 — p$) 2 
(see text for details). This generalization allows a qualitative 
description of the system for ps < 1. 

The rationale for the EW approximation (restricted to 
ps = 1) is simple 38|: in an excitable tree, activity is 
decomposed in forward- and backward-propagating ex- 
citable waves. We separate (for g > 0) the active state 
(1) into three different active states: (I A), (IP), and 
(1C), as represented in Fig. [5] P/(lvl) (black) corre- 
sponds to the density of active sites which received an 
external input and thus generates excitable-wave prop- 
agating both forwards and backwards. Pf(lP) (green) 
corresponds to the density of active sites which received 
only forward-propagating input. Finally, P t 9 (lC) (yel- 
low) corresponds to the density of active sites which re- 
ceived only backward-propagating input. 



For a correct description of the general case, ps < 1, 
we need to introduce loops into the topology as repre- 
sented by the purple arrows of Fig. O First, we assume 
that with probability 1 — p$ an active site remains ac- 
tive for the next time step. Second, among those sites 
that remained active we consider that a fraction 1 — ps 
of them can jump to P^ +1 {1A), and therefore contribute 
with both forward- and backward-propagating excitable- 
waves. This last step mimics the actions of the returning 
activity, since state {1A) is a required intermediate step 
to fulfill a back and forth movement of the excitable- wave 
under the GEW mean-field approximation. 

Following these ideas, and applying the usual mean- 
field approximations [38[ , one can write the equations for 
the g > layers as 

P/ +1 (L4) = Pf(0)Ai + (l-^) 

.{Pf(lA) + (l-tf) 

•[P/(15) + Pf(lC)]}, 
P/ +1 (lP) = P/(0)(l-A^)A^(i) 

+pl(l-pl)P t 9 (lB), 

P/ +1 (1C) = Pf(0)(l - A° A ) [1 - A|(t)] A 9 c (t) 

+pl(l-p°)Pi(lC), 

where the excitation probabilities are given by 
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Equations (|12[) and (|13j) remain unchanged, with 
P/(l) = Pf(lA) + P t a (lB) + P t 9 {lC). The dynamics 
of the most distal layer g = G is obtained by fixing 
Ag(t) = 0. The apical element (g = 0) has a simpler dy- 
namics since it does not receive backpropagating waves, 
so its activity is simply governed by 

P° +1 (l) = P°(0)A°(i) + (1 -p§)if (1) , (22) 



+Pr\lB)]}\ 



k g c {t)=p PX pr\iA)+pr\ic) 



9" 1 , 



with 



fc+i 



A°(t) = 1 - (l-p h ) {1-px [P?(1A) + P?0.B)]} 

Taking into account the normalization conditions, the 
dimensionality of the map resulting from the GEW ap- 
proximation is the same as for the EW approxima- 
tion [38(1 : 4(G — 1) + 5. Numerical solutions of this map 
leads to the solid curves in Fig. [3] and |H 
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